clear; clc;
% 固定 f = 0.7，扫描幅值 A，绘制 bifurcation diagram 和 LE vs A

% 参数设定
a = 0.7; 
b = 0.8; 
c = 0.1; 
alpha = 0.25; 
beta = 0.01; 
delta = 0.1;
% 初始状态 [x; y; z]
y0 = [0.2; 0; 0.01];

T_total = 100;   % 总模拟时间
T_trans = 20;    % 舍去瞬态时间
dt = 0.01;       % 用于 LE 估计的重新归一化间隔
h = 0.01;        % RK4 积分步长

% 固定 f=0.7，扫描幅值 A
f_fixed = 0.7;
A_range = linspace(0, 1, 100);  % 幅值扫描范围

% 用于存储分岔图数据 (A, x_samples)
bifurcationA = [];
% 用于存储 LE 值
LE_A = zeros(size(A_range));

for i = 1:length(A_range)
    A_val = A_range(i);
    params = struct('a', a, 'b', b, 'c', c, 'alpha', alpha, ...
                    'beta', beta, 'delta', delta, 'A', A_val, 'f', f_fixed);
                
    % 数值积分，使用自写的 RK4 积分器
    tspan = [0, T_total];
    [t, Y] = rk4(@(t,y) hybridNeuron(t, y, params), tspan, y0, h);
    
    % 舍去瞬态数据
    idx = t > T_trans;
    t_post = t(idx);
    Y_post = Y(idx,:);
    
    % 根据外部激励周期采样：每个周期采集 x 的峰值
    T_period = 2*pi / f_fixed; 
    x_samples = [];
    for k = ceil(T_trans/T_period):floor(T_total/T_period)
        t_sample = k * T_period;
        [~, idx_sample] = min(abs(t - t_sample));
        x_samples(end+1) = Y(idx_sample,1);
    end
    
    % 将当前 A 值对应的分岔数据 [A, x_sample] 累加起来
    bifurcationA = [bifurcationA; repmat(A_val, length(x_samples), 1), x_samples'];
    
    % 计算 LE（使用自写的 RK4 积分器）
    LE_A(i) = LE(params, y0, T_total, dt, h);
end

% 绘图
subplot(1,2,1);
plot(bifurcationA(:,1), bifurcationA(:,2), 'r.', 'MarkerSize', 5);
xlabel('A');
ylabel('x_{peak}');
title('Bifurcation diagram (f=0.7)');

subplot(1,2,2);
plot(A_range, LE_A, 'g-', 'LineWidth', 1.5); hold on;
yline(0, 'r--');  % 绘制 LE=0 的参考线
xlabel('A');
ylabel('LE');
title('LE vs A (f=0.7)');
